#!/usr/bin/perl -w
use strict;
use Getopt::Long;

my @files = qw(CJB4 CJB5-2 CJB7 JEL289 JEL627  
	       SRS812 TST75);
my @genomes = qw(JEL423 JAM81);

my $basedir = '/srv/projects/stajichlab/projects/bd_popgen/analysis';
my $datadir = '../Cornell_Run_1';
GetOptions(
	   'd|dir:s' => \$basedir,
	   'a|data:s' => \$datadir);

for my $file ( @files ) {
    for my $genome ( @genomes ) {
	open(OUT, ">map_$file.$genome.sh");
	print OUT "#!/bin/bash\n";
	print OUT "GENOME=$genome\n";
	print OUT "STRAIN=$file\n";
	print OUT "cd $basedir\n";
	print OUT "bwa aln -t 4 \$GENOME $datadir/\$STRAIN.fq > \$STRAIN.\$GENOME.sai 2> /dev/null\n";
	print OUT "bwa samse \$GENOME \$STRAIN.\$GENOME.sai $datadir/\$STRAIN.fq > \$STRAIN.\$GENOME.sam 2> /dev/null\n";
	print OUT "samtools import \$GENOME.fa.fai \$STRAIN.\$GENOME.sam \$STRAIN.\$GENOME.unsrt.bam 2> /dev/null\n";
	print OUT "samtools sort \$STRAIN.\$GENOME.unsrt.bam  \$STRAIN.\$GENOME 2> /dev/null\n";
	print OUT "samtools index \$STRAIN.\$GENOME.bam 2> /dev/null\n";
	print OUT "samtools pileup -cv -f \$GENOME.fa \$STRAIN.\$GENOME.bam > \$STRAIN.\$GENOME.pileup 2> /dev/null\n";
	print OUT "samtools.pl varFilter \$STRAIN.\$GENOME.pileup | awk \'\$6>=20\' > \$STRAIN.\$GENOME.pileup_filtered 2> /dev/null\n";
    }
}
